using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using NUnit.Framework;

namespace dnAnalytics.Tests.LinearAlgebra.Decomposition
{
    public abstract class LUTests
    {
        private const string matrixFile = "./TestData/Matrices/random_real_general_array_100.mtx";
        private const string singularFile = "./TestData/Matrices/gear_integer_general_coordinate_100.mtx";
        private LU lu;
        private Matrix matrix;
        private Matrix singluarMatrix;
        private LU singularLu;
        protected abstract Matrix GetMatrix(string file);
        protected abstract Matrix GetMatrix(int size);
        protected abstract Matrix GetIdentityMatrix(int size);
        protected abstract Vector GetVector(int size);

        [TestFixtureSetUp]
        public void SetUp()
        {
            matrix = GetMatrix(matrixFile);
            lu = new LU(matrix);
            singluarMatrix = GetMatrix(singularFile);
            singularLu = new LU(singluarMatrix);
        }

        [Test]
        [ExpectedException(typeof (MatrixNotSquareException))]
        public void LUConstructorNotSquare()
        {
            new LU(new DenseMatrix(3, 2));
        }

        [Test]
        public void Singular()
        {
            Assert.IsTrue(singularLu.IsSingular());
        }

        [Test]
        public void NotSingular()
        {
            Assert.IsFalse(lu.IsSingular());
        }

        [Test]
        public void Factors()
        {
            int[] pivots = lu.Pivots();

            Matrix perm = new DenseMatrix(pivots.Length);
            for (int i = 0; i < pivots.Length; i++)
            {
                perm[i, i] = 1;
            }

            //there has to be a quicker way to generate the perm matrix?
            double temp;
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    for (int j = 0; j < pivots.Length; j++)
                    {
                        temp = perm[i, j];
                        perm[i, j] = perm[pivots[i], j];
                        perm[pivots[i], j] = temp;
                    }
                }
            }
            Matrix lower = lu.LowerFactor();
            Matrix upper = lu.UpperFactor();
            Matrix result = lower*upper;

            Matrix expected = perm*matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    TestHelper.TestRelativeError(expected[i, j], result[i, j], Constants.AcceptableError);

                }
            }

            pivots = singularLu.Pivots();
            perm = new DenseMatrix(pivots.Length);
            for (int i = 0; i < pivots.Length; i++)
            {
                perm[i, i] = 1;
            }

            //there has to be a quicker way to generate the perm matrix?
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    for (int j = 0; j < pivots.Length; j++)
                    {
                        temp = perm[i, j];
                        perm[i, j] = perm[pivots[i], j];
                        perm[pivots[i], j] = temp;
                    }
                }
            }

            lower = singularLu.LowerFactor();
            upper = singularLu.UpperFactor();
            result = lower*upper;

            expected = perm*singluarMatrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    TestHelper.TestRelativeError(expected[i, j], result[i, j], Constants.AcceptableError);
                }
            }
        }


        [Test]
        public void Determinant()
        {
            Assert.AreEqual(0, singularLu.Determinant());
            double expected = -5.006656015927541e248;
            double actual = lu.Determinant();
            double error = System.Math.Abs((actual - expected)/expected);
            Assert.IsTrue(error < Constants.AcceptableError);
        }

        [Test]
        public void Inverse()
        {
            Matrix inverse = lu.Inverse();
            Matrix result = new DenseMatrix(inverse.Rows);
            matrix.Multiply(inverse, result);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.AcceptableError);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof (SingularMatrixException))]
        public void InverseSingular()
        {
            singularLu.Inverse();
        }
    }
}